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I. INTRODUCTION 

There is a wide range of particle systems from compu- 
tational science problems that can be efficiently and accu- 
rately simulated using asynchronous event-driven (AED) 
algorithms. Two prominent examples are molecular dy- 
namics (MD) [T] for systems of hard particles such as 
disordered granular packings [2], polymers [3J, colloids 
[4], particle-laden flows [5], and others, as well as ki- 
netic Monte Carlo (KMC) jS] for diffusion-limited reac- 
tions [7], epitaxial growth |5], quantum systems [S], bio- 
chemical reaction networks [ID], self-assembly [TT], and 
many others. As of yet unexplored are multi-scale and 
multi-physics algorithms such as event-driven dislocation 
dynamics, phase separation/precipitation, combined flow 
and diffusion with (bio)chemical reactions, etc. 

In this work we will focus on a class of particle-based 
problems that are very common in computational ma- 
terials science and are well-suited for AED simulation. 
Specifically, we will focus on the simulation of large sys- 
tems of mobile particles interacting with short-range pair- 
wise (two-body) potentials (forces). Our goal will be 
to reveal the common building blocks of these simula- 
tions (e.g., event queues, neighbor searches), but also to 
highlight the components that are problem specific (e.g., 
event prediction and processing). We will present these 
components in some detail for three specific examples: 
event-driven molecular dynamics (EDMD), first-passage 
kinetic Monte Carlo (FPKMC), and stochastic EDMD 
(SEDMD). Through the discussion of these examples we 
will demonstrate the undeniable advantages of AED algo- 
rithms, but we will also reveal the difficulties with using 
AED algorithms for realistic models. 



A. Background 

We consider the simulation of the time evolution of 
a collection of N interacting particles in (i-dimensions, 
starting from some initial condition. At any point in 
time, the system Q — (Q, B) is characterized by the con- 



figuration Q = (qi, . . . , qjv), containing at least the cen- 
troid positions r 4 for every particle i, and the additional 
global information B, which may involve variables such 
as boundary conditions or external fields. The number of 
particles N may itself vary with time. For each particle i 
we may consider an arbitrary number of attributes in 
addition to the position of the centroid r<, q^ = (r^a*), 
for example, a^ may also contain the linear and/or angu- 
lar velocity, the orientation and/or the chemical species 
(shape, charge, mass, internal composition) of particle i. 
Typically a^ will at least contain an integer that identi- 
fies its species 1 < s, < N Sl and some information will be 
shared among all particles belonging to the same species 
(ex., the charge or mass). In particular, the symmetric 
interaction table T a p stores N S (N S + l)/2 logical (true or 
false) entries that specify whether species a and (3 inter- 
act or not. 

Two particles i and j are overlapping only if a certain 
(generalized) distance between them (fy(qi,qj) > is 
less than some cutoff distance or diameter Dij. Over- 
lapping particles react with each other in an application 
specific manner. Typically the type of reaction and fiy 
depend only on the species of the two particles, but there 
may also be dependencies on time or some other exter- 
nal field parameters. For example, for (additive) hard 
spheres = ||rj — Tj\\ L is the L 2 norm of = — r\,-, 
and = (Di + Dj)/2 and the type of reaction is (hard- 
core) repulsion. For a non-interacting pair of particles 
one may set < 0. Particles may also overlap with 
boundaries of the simulation domain, such as hard walls 
or reactive surfaces, however, typically the majority of in- 
teractions are among particles. Our assumption of short- 
range interactions implies that all Dij 's are much smaller 
than the system size. If there are long range interac- 
tions present (e.g., electrostatic forces in plasmas), it may 
sometimes be possible to split them into a short-range 
part and a long-range part and treat the long-range part 
separately (e.g., multipole or cell-based FFT methods) 
as part of the "boundary" handling. 

We will assume that the time evolution (motion) of the 
system is mostly smooth with the exception of certain dis- 
crete events, which lead to discontinuous changes in the 
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configuration of a certain collection of particles. Events 
may involve a single particle i, such as the change of the 
internal state of the particle (e.g., decay reactions, spin 
flips, sudden changes in the particle velocity). Events 
may also involve pairs of particles, for example, the col- 
lision (exchange of momentum) between two particles i 
and j, or a chemical reaction between two overlapping 
particles i and j leading to the creation and/or destruc- 
tion of particles. For now we will ignore events involving 
more than two particles. Some events may also involve 
global variables in B and thus implicitly affect all of the 
particles. Here we will assume events are instantaneous, 
i.e., they have no duration. Events that have a duration 
(e.g., particles overlapping for a certain time interval) 
can be treated as a pair of events, one for the start of the 
event and one for the end of the event. In a more general 
framework one may consider all events as having a finite 
duration, as in process-oriented simulations |12j (in such 
a framework collisions would be a special degenerate case 
of a more general "overlap" event). 

We will refer to simulations as event-driven (also called 
discrete event simulations |13j ) if the state of the sys- 
tem is evolved discontinuously in time from one event 
to another, predicting the time of the next event when- 
ever an event is processed. This is in contrast with the 
more common time-driven simulations, where the state 
changes continuously in time and is evolved over a se- 
quence of small time steps At, discovering and process- 
ing events at the end of the time step. Time-driven al- 
gorithms are inaccurate when the time step is large, and 
they become asymptotically accurate as At — > 0. There- 
fore, the time step must be smaller than the estimated 
fastest time scales in the problem. This leads to large in- 
efficiencies when there are multiple dynamically-changing 
time scales. On the other hand, event-driven algorithms 
automatically adjust the time step. 

We will focus on asynchronous event-driven (AED) al- 
gorithms. In asynchronous algorithms, there is a global 
simulation time t, typically the time when the last pro- 
cessed event occurred, and each particle is at a different 
point in time ti < t, typically the last time it participated 
in an event. This is to be contrasted to synchronous 
event-driven algorithms, where all of the particles are at 
the same time t. One of the most important examples 
of a synchronous event-driven algorithm in materials sci- 
ence is the n-fold (BKL) algorithm for performing kinetic 
(dynamic) Monte Carlo simulations [5]. The exactness 
and efficiency of the n-fold algorithm hinge on the fact 
that the state of the system does not change in-between 
events, as is common in lattice models where the posi- 
tions of the particles are discrete. For example, atoms 
may stay in the immediate vicinity of the crystal lattice 
sites and only sometimes hop to nearby sites. In the types 
of problems that we will consider, the positions of the 
particles will be continuous and continuously changing 
even in-between events. Therefore, synchronous KMC 
algorithms will be approximate and can be viewed as the 
equivalent of time stepping, where the time step At is 



not constant. It is important to point out that even 
in cases where the evolution of the system consists en- 
tirely of discrete jumps (e.g., Markov chain transitions), 
asynchronous algorithms may be more efficient than syn- 
chronous ones [IT]. It is possible to combine the two 
algorithms by using the more general asynchronous algo- 
rithm at the top level but treating a subset of the par- 
ticles synchronously, as if they are a super-particle with 
complex internal structure. 



B. Model Examples 

Atomistic or molecular-level modeling is one of the 
foundations of computational materials science. The two 
most popular types of algorithms used in the simulation 
of materials are molecular dynamics (MD) and Monte 
Carlo (MC) algorithms. Monte Carlo algorithms are of- 
ten used to study static equilibrium properties of systems, 
however, here we focus on dynamic or kinetic Monte 
Carlo, where the time evolution of a system is modeled 
just like in MD. For our purposes, the only important 
difference between the two is that MD is a deterministic 
algorithm, in which deterministic equations of motions 
are solved, and Monte Carlo is a stochastic procedure in 
which sample paths from an ensemble of weighted paths 
are generated. In both cases one typically averages over 
multiple trajectories, starting with different initial condi- 
tions and/or using a different random number seed. From 
the perspective of AED algorithms, this means that ran- 
dom number generators (RNGs) are involved in the de- 
termination of the time certain events occur as well as in 
the actual processing of those events. 

The very first molecular dynamics (MD) calculations 
simulated a system of hard disks and hard spheres and 
used an AED algorithm [JJ. Event-driven MD (EDMD) 
algorithms for hard particles are discussed in considerable 
detail in Ref. [2] and elsewhere, here we only present the 
essential components. The hard-sphere system is a col- 
lection of non-overlapping spheres (or disks), contained 
within a bounded region, each moving with a certain ve- 
locity Vj = r'j. Pairs of spheres collide and the collid- 
ing particles bounce off elastically, preserving both lin- 
ear momentum and energy. Many generalizations can be 
considered, for example, the spheres may be growing in 
size and/or the particles may be nonspherical [14], the 
collisions may not be perfectly elastic [TS], some of the 
particles may be tethered to each other to form structures 
such as polymer chains [3], etc. The general features are 
that particles move ballistically along simple determinis- 
tic paths (such as straight lines) in-between binary col- 
lisions. The primary type of event are binary collisions, 
which have no duration and involve deterministic changes 
of the velocities of the colliding particles. The ballistic 
motion of the particles is described by Newton's equa- 
tions of motion (i.e., deterministic ODEs), mv; = Fj, 
where F is an external forcing (e.g., gravity). 

Direct simulation Monte Carlo (DSMC) [T5] is an MC 
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algorithm that tries to mimic MD for fluids. We will con- 
sider DSMC as a fast alternative to MD, even though it 
can also be viewed as a particle-based MC method for 
solving the Boltzmann equation in dilute fluids. From 
our perspective, DSMC is an approximate variant of MD 
in which the particle collisions are not processed exactly, 
rather, particle collisions are stochastic and (attempt to) 
follow the same probability distributions as would have 
exact MD. Specifically, nearby particles are randomly 
chosen to undergo stochastic collisions and exchange mo- 
mentum and energy, thus leading to local conservation 
laws and hydrodynamic behavior. DSMC is applicable 
in cases when the structure of the fluid and the detailed 
motions of all of the particles do not matter, as is the 
case with solvent molecules (e.g., water) in fluid flow 
problems or large-scale granular flows |15j . Tradition- 
ally DSMC has been implemented using a time-driven 
approach, in which at each time step particles are first 
propagated in a ballistic (convective) fashion, and then 
a certain number of stochastic particle collisions among 
nearby particles are processed. Here we describe a novel 
AED algorithm for DSMC, and demonstrate how it can 
be integrated with EDMD in order to replace the expen- 
sive MD with cheaper DSMC for some of the particle 
species (e.g., solvent molecules). We term the resulting 
algorithm Stochastic EDMD (SEDMD). 



The motion of the particles in-between events is not 
always deterministic. In particular, an important class 
of problems concerns diffusing particles, that is, particles 
whose velocity changes randomly very frequently (i..e, 
they make many small steps in random directions). The 
motion of the particles is probabilistic, in the sense that 
the probability c(r, t + At) of finding a particle at a given 
position r at a certain time t + At, assuming it started at 
the (space and time) origin, is the solution to the time- 
dependent diffusion equation dtc = DV 2 c (a determinis- 
tic partial differential equation), where D is the particle 
diffusion coefficient (generally a tensor) . A variety of re- 
actions may occur when a pair of particles collides, for ex- 
ample, particles may repel each other (colloids [4 ), they 
may stick or begin merging together (paint suspensions) , 
or they may undergo a chemical reaction that consumes 
the reacting particles and produces zero, one, two, or 
possibly more new particles (a wide range of diffusion- 
limited reactions in materials [Tj [T7] and biological sys- 
tems p~7]). Several approximate event-driven KMC algo- 
rithms have been used in the past for this problem |17j . 
Here we will describe a recently-developed First-Passage 
Kinetic Monte Carlo (FPKMC) exact AED algorithm for 
simulating a collection of diffusing hard particles [7]. It 
is worth pointing out that one may also consider par- 
ticles whose trajectories are a combination of ballistic 
and diffusive motion, that is, motion that is described by 
Langevin's equations (or other stochastic ODEs or even 
PDEs). In that sense, we will see that both the MD and 
MC algorithms share many common features. 



II. A GENERAL AED PARTICLE ALGORITHM 

In this work we focus on systems where particles only 
interact with nearby particles. We will formalize this by 
defining a geometric hierarchy of regions around a given 
particle. These particle proximity hierarchies are at the 
core of geometry-specific (GS) aspects of AED simulation, 
which can be reused for different application-specific (AS) 
rules for moving and interacting the particles. We will 
assign a hard core Ci to each particle such that a parti- 
cle may overlap with another particle only if their cores 
overlap. For (additive) hard spheres, the core is nothing 
more than the particle itself. Next, we protect particle 
i against other particles by enclosing it inside a protec- 
tive region Vi, Ci C "Pi, that is typically disjoint from 
the majority of other protective regions. Finally, we as- 
sume that every protective region i is contained within a 
neighborhood region Mi, Vi C Mi. The set of neighbors of 
i consists of the particles j whose neighborhood regions 
intersect Mi, Mi fl Mj ^ 0, and which are of a species 
interacting with the species of particle i, i.e., T SiS . = T. 

We will assume that when a particle does not interact 
with other particles we can easily follow its time evo- 
lution (motion), that is, given the current configuration 
Qi(t), we can probabilistically determine the position at 
a later time qj(f + At). This is a single-particle problem 
and can typically be solved analytically. For example, in 
MD the particle trajectory is a unique (i.e., deterministic) 
straight path, r,(t+At) = rj(i) + v^Ai, while in diffusion 
problems it is the solution to a (stochastic) Langevin or 
diffusion equation. Event-driven algorithms are efficient 
because they use such analytic solutions to quickly prop- 
agate particles over potentially large time steps as long 
as they are far enough from other (interacting) particles. 
We will also assume that one can solve two-body prob- 
lems for the case when two particles are isolated from 
other particles but may interact with each other. These 
two-body problems are typically much more difficult to 
solve (quasi) analytically. Specific examples will be given 
later. 



A. The Event Loop 

An AED algorithm consists of processing a sequence 
of time-ordered events. Each particle i must store some 
basic information needed to predict and process events 
associated with it. The particle time ti specifies the last 
time the configuration of particle i was updated, tj < t, 
where t is the current simulation time. Some particles 
may be time-driven and thus not have their own event 
prediction. The rest of the particles are event-driven and 
each such particle stores a prediction for its impending 
event (t e ,p,v), specified via the predicted time of oc- 
currence t e (a floating point number), the event partner 
p (an integer), and the event qualifier (type of event) 
v (also an integer). The partner p could be some pre- 
specified invalid value to identify time-driven particles. 



4 



Note that the event schedules must be kept symmetric 
at all times, that is, if particle i has an impending event 
with j, then particle j must have an impending event 
with i. A particle may store multiple event predictions, 
in order to avoid re-predicting events if the impending 
event is invalidated, however, we will not explicitly han- 
dle this possibility due to the complications it introduces. 

The exact interpretation of p and v, for a given particle 
i, is application- and geometry-specific. Some common 
types of events can be pre-specified by reserving certain 
values of the event partner p, for example, we have used 
the following list for the set of models presented here: 

p = An update of the event prediction for i, not re- 
quiring an update of q^. The value v = means 
that q.; has not changed since the last prediction for 
i (thus allowing stored information from previous 
predictions to be reused if needed), v = 1 means 
that an event occurred which did not alter the ge- 
ometry (for example, the position of i is the same 
but its velocity changed), while v = — 1 means that 
this particle was just inserted into the system and 
a geometry update is necessary as well. 

p = i > A single-particle event that requires an update 
of qi . The special value v = denotes a simple time 
advance of i without any additional event process- 
ing, v < denotes an event that does not change 
the geometry (for example, only the velocity of a 
particle changes), and v > is used for additional 
AS events that may also change the geometry (e.g., 
particle decay). 

1 < p < N and p i An unprocessed binary reaction 
between particles i and j = p, with additional AS 
information about the type of reaction stored in 
v, for example, elastic collision, a certain chemical 
reaction, etc. 

p = oo A "boundary" event requiring the update of the 
particle geometry. If v = then only the protective 
region Vi needs to be updated, if v = — 1 then 
the neighborhood Mi needs to be updated (collision 
with a virtual boundary) , v < — 1 denotes collisions 
with pre-specified boundaries (such as hard walls), 
and v > specify AS boundary events (such as 
collisions with reactive surfaces). 

p = — oo Denotes an invalid event, meaning that this par- 
ticle is not in the event queue and is handled sepa- 
rately, for example, it is time-driven. 

It is important to point out that we are not suggesting 
that an actual implementation needs to use integers to 
identify different types of events. In an object-oriented 
framework events may be represented as objects that in- 
herit from a base event class and have methods to handle 
them, with the base implementation providing handlers 
for certain pre-defined (single, pair, and boundary) types 



of events. We do not discuss here the possible ways to or- 
ganize an inheritance hierarchy of classes for AED simu- 
lations, since such a hierarchy involves multiple complex 
components, notably a module for handling boundary 
conditions in static and dynamic environments, a module 
for handling static and dynamic particle geometry (over- 
lap, neighborhoods, neighbor searches, etc.), an event- 
dispatcher, a visualization module, application-specific 
modules for event scheduling and handling, etc. 

Algorithm [l] represents the main event loop in the 
AED algorithm, which processes events one after the 
other in the order they occur and advances the global 
time t accordingly. It uses a collection of other auxiliary 
geometry-specific (GS) or application-specific (AS) steps, 
as marked in the algorithm outline. Specific examples of 
various GS and AS steps are given in the next section. 



Algorithm 1: Process the next event in the event 
queue. 



1. Find (query) the top of the event queue (usually a heap) 
to find the next particle i to have an event with p at t e . 
Note that steps marked as (AS+C) below may reorder 
the queue and/or cycle back to this step. 

2. Find the next "external" event to happen at time t ex , 
possibly using an additional event queue (AS). 

3. If t ex < t e then process the external event (AS) and 
cycle back to step [l] 

4. Remove i from the event queue and advance the global 
simulation time t <— t e . 

5. If p = and v = — 1 then build a new Mi and then check 
if i overlaps with any of its new neighbors (GS). If it 
does, process the associated reactions (AS+C), other- 
wise build a new Vi as in step |8a| 

6. Else if p = i then update the configuration of particle 
i to time t using a single-particle propagator (AS), and 
set ij <— t. If v 7^ then process the single-particle 
event (AS+C). If v > then search for overlaps as in 
step [5] 

7. Else if 1 < p < N then update the configuration of 
particles i and j = p using a two-particle propagator 
(AS), set ti*—t and tj <— t, and then process the binary 
reaction between i and j (AS+C). This may involve 
inserting particle j back into the queue with t e = t, 
p = 0,v = 0. 

8. Else if p = oo, then update and ti as in step [6] 
If v > then process the boundary event (AS+C), 
otherwise 

(a) If v — then update Vi (AS+GS), typically in- 
volving an iteration over the neighbors of i. 

(b) Else if v = — 1 then update Mi and identify the 
new neighbors of particle i (GS). 

(c) Else if v < —1 process the geometry-induced 
boundary event (GS+AS). 

9. Predict a new t e , p, and v for particle i by finding the 
minimal time among the possible events listed below. 
Each successive search needs to only extend up to the 
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current minimum event time, and may return an incom- 
plete prediction t e > t, p — 0, v = 0, where t e provides 
a lower bound on the actual event time. 

(a) When particle i leaves Vi or A/i (AS). 

(b) When particle i undergoes a single-particle event 
(AS). 

(c) When particle i first reacts with a neighbor j 
(AS), as found by searching over all neighbors j 
whose protective region Vj intersects Vi (GS). If a 
particle j gives the current minimum event time, 
remove it from the event queue. If such a parti- 
cle j has an event partner that is another parti- 
cle (third party) k =fc i, update the positions of 
j and k using the two-particle propagator as in 
step [7j invalidate fc's event prediction by setting 
its t e <— t, p <— 0, v <— 0, and update its position 
in the event queue (alternatively, one may use lazy 
invalidation strategies) . 

10. Insert particle i back into the event heap with key t e 
and go back to step[l] 



1. Non-Particle Events 

In a variety of applications the majority of events can 
be associated with a specific particle, and one can sched- 
ule one event per particle in the event queue. However, 
sometimes there may be events that are associated with 
a (possible large) group of particles, or events that are 
not specifically associated with a particle. We consider 
these non-particle events as application-specific "exter- 
nal" events in Algorithm [T] 

An important example of such an event are time step 
events. Namely, some group of particles may not be 
propagated asynchronously using the event queue, in- 
stead, the particles in the group may be updated syn- 
chronously, for example, in regular time intervals At. 
There may in fact be multiple such groups each with 
their own timestep, for example, each species might have 
its own time step. Alternatively, all or some of the parti- 
cles may be updated in a time-stepped manner and addi- 
tional asynchronous events may be processed in-between 
the time step events. An example is molecular dynamics 
in which time-driven handling is combined with event- 
driven handling for the hard-core collisions [5J |TH . The 
time step events should also be ordered in time and the 
next one chosen as the external event in Step [2] in Al- 
gorithm [T] A separate priority queue may be used for 
ordering the external events. In general, there may be 
an event queue of events associated with particles, with 
cells, with species, etc. These may be separate queues 
that are joined at the top or they can be merged into a 
single heap. 



B. Near-Neighbor Search 

All large-scale particle-based algorithms use various ge- 
ometric techniques to make the number of neighbors of a 
given particle O(l) instead of O(N). Reference [5] pro- 
vides extensive details (and illustrations) of these tech- 
niques for hard spheres and ellipsoids; here we summarize 
only the essential components. 

1. Linked List Cell (LLC) Method 

The most basic technique is the so-called linked list cell 
(LLC) method. The simulation domain, typically an or- 
thogonal box, is partitioned into N c cells, typically cubes. 
Each particle i stores the cell q to which its centroid be- 
longs, and each cell c stores a list C c of all the particles it 
contains (usually we also store the total number of par- 
ticles in the cell). Given a particle and a range of inter- 
action, the lists of potential neighbors is determined by 
scanning through the neighboring cells. For maximal ef- 
ficiency the cell should be larger than the largest range of 
interaction so that only the nearest-neighbor cells need 
to be searched. There are more sophisticated neighbor 
search methods developed in the computational geom- 
etry community, such as using (colored) quad/oct-trees, 
however, we are not aware of their use in AED implemen- 
tations, likely because of the implementation complexity. 
This is an important subject for future research. 

It is important to point out that in certain applications 
the cells themselves play a crucial role in the algorithm, 
typically as a means to provide mesoscopic averages of 
physical variables (averaged over the particles in a given 
cell) used to switch from a particle-based model to a 
continuum description. For example, in PIC (particle- 
in-cell) algorithms for plasma simulation, the cells are 
used to solve for background electric fields using FFT 
transforms [19 . In DSMC, the algorithm stochastically 
collides pairs of particles that are in the same cell. In 
some applications, events may be associated with the 
cells themselves, instead of or in addition to the events 
associated with particles [20]. Usually the same cells are 
used for both neighbor searches and averaging for sim- 
plicity, however, this may not be the optimal choice in 
terms of efficiency. 

For a method that only uses the LLC method for neigh- 
bor searches, the neighborhood region A/", is composed of 
the (typically 3 d , where d is the dimensionality) cells that 
neighbor a, including a itself. The protection region Vi 
may be a simple geometric region like a sphere inscribed 
in A/j (sphere of diameter smaller than the cell size), it 
may be that Vi = Cj + C{ , or maybe Vi =N%. 

2. Near-Neighbor List (NNL ) Method 

Another neighbor search method is the near-neighbor 
list (NNL) method, which is described for hard particles 
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in Ref. [SJ. The idea is to use as Mi a region that (when 
it is created) is just an enlargement of the particle by 
a certain scaling factor /i > 1. When Mi is created the 
method also creates a (linked) list of all the neighbor- 
hoods that intersect it, to form NNL(i) (hard walls or 
other boundaries may also be near neighbors). This list 
of (potential) interactions can then be reused until the 
particle core Ci protrudes outside of Mi. This reuse leads 
to great savings in situations where particles are fairly 
localized. 

Note that the LLC method is still used in order to 
create Mi and NNL(i) even if NNLs are used, in order 
to keep the maximal cost of pairwise searches at O(N) 
instead of 0(N 2 ). In some situations (such as mixed 
MD/DSMC simulations as we describe later) one may 
use NNLs only for a subset of the particles and use the 
more traditional LLCs for others. In this case one can 
use Vi = Mi n (ci + Ci) in order to ensure that both the 
NNLs (for those particles that have them) and the LLCs 
(for all particles) are valid neighbor search methods. 



3. Cell Bitmasks 

Efficient handling of spatial information is an essential 
component of realistic AED algorithms. For example, 
further improvements to the basic LLC method may be 
required for certain applications in order to maximize the 
efficiency of neighbor searches. The handling of boundary 
conditions or domain-decomposition is an application- 
specific component that is in some sense disjoint from 
the basic AED framework presented here, however, it is 
often very important in practice. In this section, which 
may be skipped at first reading, we describe an enhance- 
ment to the basic LLC method that we have found very 
useful in handling spatial information. 

In our implementation, in addition to the list of par- 
ticles C c , each cell c stores a bitmask M. c consisting of 
•^bits > ^ s kits (bitfields), where N s is the number of 
species. These bits may be one (set) or zero (not set) to 
indicate certain properties of the cell, specifically, what 
species of particles the cell contains, whether the cell is 
event or time driven, and to specify boundary conditions. 
Bit 7 in the bitmask A4 C is set if the cell c contained a 
particle of species 7 in the near past and may still con- 
tain a particle of that species. The bit is set whenever a 
particle of species 7 is added to the cell, and it should be 
cleared periodically if the cell no longer contains parti- 
cles of that species. When performing a neighbor search 
for a particle i, cells not containing particles of species 
that interact with species Si are easily found (by OR'ing 
the cell masks with the Sj'th row of T SiSj ) and are sim- 
ply skipped. This can significantly speed up the neighbor 
searches in cases where not all particles interact with all 
other particles. 

For the purposes of a combined event-driven with time- 
driven algorithm one may also need to distinguish those 
cells where particles are treated using a fully event-driven 



(ED) scheme. We use one of the bits in the bitmasks, bit 
lEDi to mark event-driven (ED) cells, and the choice 
of such cells is in general application specific. For ex- 
ample, for diffusion problems, cells with a high density 
may be treated more efficiently using time-driven (small 
hopping) techniques while areas of low density may be 
treated more efficiently using the asynchronous event- 
driven algorithm. In our combined MD with DSMC algo- 
rithm cells that contain or that are near non-DSMC par- 
ticles are event-driven while those containing only DSMC 
particles are treated as time-driven cells. Note that all of 
the cell bitmasks should be reset and then re-built (i.e., 
refreshed) periodically 

The cell masks can also be used to specify partitionings 
of the simulation domain. This is very useful in specify- 
ing more general boundary conditions in situations when 
the event-driven simulation is embedded inside a larger 
domain that is not simulated explicitly. For example, a 
molecular dynamics simulation may be embedded in a 
multiscale solver where the surrounding space is treated 
using a continuum method (finite element or finite vol- 
ume, for example) coupled to the particle region through 
an intermediary boundary layer. Similar considerations 
apply in parallelization via domain decomposition, where 
the simulation domain simulated by a single processing 
element (PE) or logical process (LP) is embedded inside a 
larger region where other domains are simulated by other 
PEs/LPs. 

We classify the cells as being interior, boundary, and 
external cells (a specific illustrative example is given in 
Fig. [lj. Our implementation uses bits in the cell bit- 
masks to mark a cell as being boundary (bit 7s), or 
external (bit 7^), the rest are interior. Interior cells 
are those for which complete boundary conditions are 
specified and that cannot be directly affected by events 
occurring outside of the simulation domain (the interior 
cells are divided into event-driven and time-driven as dis- 
cussed earlier) . Boundary cells surround the interior cells 
with a layer of cells of thickness w B > 1 cells (typically 
w B = 1 or w B = 2) and they represent cells that are 
affected by external events (i.e., events not simulated 
directly). External cells are non-interior cells that are 
not explicitly simulated, rather, they provide a bound- 
ary condition / padding around the interior and boundary 
cells. This layer must be at least w B cells thick, and the 
cells within a layer of w B cells around the simulation do- 
main (interior together with boundary cells) are marked 
as both external and boundary cells (e.g., these could be 
ghost cells in parallel simulation) . It is often the case that 
w B = Wg = wb- All of the remaining cells are purely 
external cells and simply ignored by the simulation. 



III. MODEL EXAMPLES 

In this section we present the handling of the various 
AS and GS steps in Algorithm [T] for three specific model 
applications. 
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A. Event-Driven Molecular Dynamics (EDMD) 

Hard-particle molecular dynamics is one of the first 
applications of AED algorithms in computational science, 
and is discussed in more detail in Ref. [2j and references 
therein. Hard-sphere MD has been used extensively in 
simulations of physical systems over the past decades and 
is also one of the few particle AED algorithm that has 
been studied in the discrete-event simulation community, 
as the billiards problem. Because of this rich history 
and extensive literature we only briefly discuss EDMD 
focusing on how it fits our general framework. 

The basic type of event in EDMD are binary colli- 
sions, which alter the momenta of two touching (and ap- 
proaching) particles, typically based elastic collision laws 
(conservation of momentum and energy) . Collisions are 
assumed to have no duration and (very unlikely) triple 
collisions are broken up into a sequence of binary ones. 
In-between collisions particles move ballistically along 
simple trajectories such as straight lines (force-free mo- 
tion) or parabolas (constant-acceleration motion). For 
hard spheres, event time predictions are based on (alge- 
braic) methods for finding the first root of a polynomial 
equation (linear, quadratic or quartic |18]L For parti- 
cle shapes that include orientational degrees of freedom, 
such as ellipsoids, numerical root finding techniques need 
to be used [SJ. 



1. AED Implementation of EDMD 

When LLCs are used, the main type of boundary event 
are cell transfers, which occur when the centroid of a par- 
ticle i collides with the boundary of the cell a . If the cell 
is at the boundary, a unit cell change occurs for periodic 
boundaries (i.e., wrapping around the torus), and hard- 
wall collisions occur for hard-wall boundaries. When 
NNLs are used, cell transfers do not have to be processed 
(i.e., one can set Vi = Mi), unless it is important to have 
accurate LLCs [as in DSMC, where V l = JVj H (c* + C t )}. 

The main AS steps in Algorithm [I] for (classical) MD 
simulations are: 

• Step [6] consists of updating the particle position, 
and possibly also velocity, ballistically. 

• Step [7] consists of updating the positions of each 
of the particles separately, as in step [6] and then 
updating their velocities taking into account the 
collisional exchange of momentum. 

• For collisions with hard walls in Step [8] the particle 
velocity is updated accordingly. Cell transfers in 
step [8c] consist of updating the LLCs by removing 
the particle from its current cell and inserting it 
into its new cell (found based on the direction of 
motion of the centroid). If the new cell is across 
a periodic boundary, the centroid is translated by 
the appropriate lattice cell vector (if NNLs are used, 



this may require updating information relating to 
periodic boundaries for each interaction). 

Step [9] is the most involved and time consuming: 



— In step 9a the time of the next cell transfer 



B. 



is predicted based on the centroid velocity. If 
NNLs are used, then the time of (virtual) col- 
lision of the core Ci with the interior wall of 
Afi is also calculated. 

In step [9c] predictions are made for binary 
collisions between particle i and each of the 
particles in neighboring cells, or between i and 
each of the particles in NNL(i). 



Stochastic Event-Driven Molecular Dynamics 
(SEDMD) 



We have recently developed a novel AED algorithm 
for simulating hydrodynamics at the molecular level that 
combines DSMC, which is a method for simulating hy- 
drodynamic transport, with event-driven molecular dy- 
namics (EDMD). The algorithm replaces some of the de- 
terministic collisions in EDMD with stochastic collisions 
and we term it Stochastic EDMD (SEDMD). The algo- 
rithm is described in more detail in Ref. [21 . 

First we describe how to transform the traditional 
time-driven DSMC algorithm [TB] into an event-driven 
algorithm, and then combine this algorithm with EDMD. 
Finally, we explain how to achieve higher efficiency by re- 
verting the DSMC component to time-driven handling. 



1. DSMC for Hydrodynamics 

The DSMC algorithm can be viewed as an approxi- 
mation to molecular dynamics in cases when the internal 
structure of the fluid, including the true equation of state, 
is not important. In particular, this is the case when sim- 
ulating a solvent in applications such as the simulation of 
large polymer chains in solution. The exact trajectories 
of the solvent particles do not really matter, and what 
really matters are the (long time and long range) hydro- 
dynamic interactions that arise because of local energy 
and momentum exchange (viscosity) and conservation 
(Navier-Stokes equations). Any method that simulates 
the correct momentum transfer localized at a sufficiently 
small scale is a good replacement for full-scale MD, and 
can lead to great computational savings when a large 
number of solvent molecules needs to be simulated. 

DSMC [16] achieves local momentum exchange and 
conservation by performing a certain number of stochas- 
tic collisions between randomly chosen pairs of particles 
that are inside the same cell. The collision rate inside a 
cell containing Nt, particles is proportional to Nl(Nl-I) 
with a pre-factor that can be based on theory or fitted to 
mimic that of the full MD simulation. For hard spheres, 
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the probability of choosing a particular pair ij is pro- 
portional to the relative velocity , and typically a re- 
jection technique (null-method technique) is used when 
choosing pairs. Specifically, the collision rate is made pro- 
portional to the maximal possible relative velocity v ™f x , 
and a randomly chosen pair ij is rejected or accepted 
with probability v^/v^ 1 ^- A pair rejection involves a 
small calculation and a random number generation and is 
thus rather inexpensive, as long as the acceptance prob- 
ability is not too small, which can typically easily be 
achieved by a judicious (but still rigorous) choice for 
u™. Alternative rules for selection of collision part- 
ners and for the choice of the post-collisional velocities 
can give thermodynamically-consistent non-ideal DSMC 
fluids [12] . 

Traditionally DSMC is performed using a time-driven 
approach: The particles are first propagated ballistically 
by a certain time step At and then sorted into cells ac- 
cordingly, and then an appropriate number of stochastic 
collisions are carried out in each cell. Furthermore, the 
existence of a finite step leads to errors in the transport 
coefficients (such as shear viscosity) of order At 2 [33], in 
addition to the errors inherent in DSMC that are of or- 
der Ax 2 [23], where Ax is the size of the cells. Therefore, 
in traditional DSMC At should be small enough so that 
particles move only a fraction of the cell size and a frac- 
tion of the mean free path during one step [53J H3]. As 
a consequence, only a fraction (10-25%) of the particles 
actually undergo a collision, and the rest of the parti- 
cles are propagated needlessly. We note, however, that 
DSMC can be extended into the dense-fluid regime in 
which case the collisional frequency is high and thus most 
particles undergo a stochastic collision at every time-step 



2. AED Implementation of DSMC 

An alternative event-driven implementation of DSMC 
explicitly predicts and process cell transfers, just as in 
EDMD algorithm. Particles positions are thus only up- 
dated when needed, and there is no time step error since 
it is easy to obtain and numerically evaluate closed-form 
expressions for the time of the cell crossings. The DSMC 
particles are represented as a species 6 for which par- 
ticles do not interact with other DSMC particles (i.e., 
Igg = J 7 ), so that the MD algorithm does not predict bi- 
nary collisions for the DSMC particles. Instead, stochas- 
tic binary collisions are added as an external Poisson 
event of the appropriate rate. One approach is to main- 
tain a global time-of-next-DSMC-collision t sc to deter- 
mine when a stochastic collision is attempted, and to 
use cell rejection to select a host cell for the collision. 
The rate of DSMC collisions is chosen according to the 
cell with maximal occupancy A r " lax , and a randomly 
chosen cell of occupancy Nl is accepted with probabil- 
ity N L (N L - l)/[iV^^(iV£ mx - 1)]. It is possible to 



avoid cell rejections altogether, using any of the multiple 
"rejection-free" techniques common to KMC simulations 
[5J. For example, the cells could be grouped in lists based 
on their occupancy and then an occupancy chosen first 
(with the appropriate weight), followed by selection of a 
cell with that particular occupancy. 

Most of the AS steps in Algorithm [T] for DSMC are 
shared with MD. The different steps are associated with 
the processing of stochastic collisions: 

• In Step [2] t ex is the time of the the next DSMC 
collision. If a rejection-free technique is used the 
cell at the top of the cell event queue is chosen, 
otherwise, a cell is selected randomly and accepted 
or rejected based on its occupancy. If the cell is 
rejected, no collisions are processed. 

• In Step [3] a pair of particles i and j is randomly se- 
lected from the previously-chosen cell, and accepted 
or rejected for collision based on the relative veloc- 
ity. If accepted, a randomly-chosen amount of mo- 
mentum and energy is exchanged among the par- 
ticles and they are moved to the top of the event 
queue with t e = t, p = 0, and v = 1. 

We have validated our event-driven DSMC algorithm by 
comparing against published results for plane Poiseuille 
flow of a rare gas obtained using traditional time-driven 
DSMC [25]. 

3. Stochastic Event-Driven Molecular Dynamics (SEDMD) 

The event-driven DSMC algorithm has few advantages 
over a time-driven approach, which are outweighed by the 
(implementation and run-time) cost of the increased algo- 
rithmic complexity. However, the AED variant of DSMC 
is very similar to EDMD, and therefore it is relatively 
simple to combine DSMC with EDMD in an event-driven 
framework. This enables the simulation of systems such 
as colloids or hard-sphere bead-chain polymers [3] in so- 
lution, where the solute particles are treated using MD, 
and the less-important solvent particles are treated ap- 
proximately using DSMC. The solvent-solute interaction 
is still treated with MD. Similar studies have already been 
carried out using time-driven MD for the solute parti- 
cles and a simplified variant of DSMC that approximates 
the solute particles as point particles and employs multi- 
particle stochastic collisions [26 . 

We have designed and implemented such a com- 
bined algorithm, which we term Stochastic Event-Driven 
Molecular Dynamics (SEDMD) [21 . The implementa- 
tion is almost identical to classical hard-sphere MD, with 
the addition of a new DSMC species 5 for which parti- 
cles do not interact with particles of the same species 
(i.e., Xss = That is, the DSMC "hard spheres" freely 
interpenetrate each other, but collide as usual with other 
species. An external stochastic collision event occurring 
as a Poisson process of the appropriate rate is used to 
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collide randomly chosen pairs of nearby DSMC particles, 
as described in the previous section. 

The algorithm is much more efficient than EDMD be- 
cause of the replacement of deterministic collisions with 
stochastic ones; however, it is still not as efficient as 
classical time-driven DSMC, especially at higher colli- 
sion rates. The two main causes of this are the over- 
head of the event queue operations in the AED variant 
of DSMC (note that the queue needs to be updated after 
every stochastic collision) , and also the cost of neighbor 
searches. Namely, for each DSMC particle the nearby 
cells need to be searched in Step [9c] to make sure they do 
not contain any solute particles (recall that cell bitmasks 
are used to efficiently implement this). In cases when 
most of the cells contain only DSMC particles, this can 
introduce significant overheads. 

This inefficiency can be corrected by combining time- 
driven with event-driven handling. Specifically, only 
those those cells that neighbor cells that contain non- 
DSMC particles are marked as event-driven (ED), the 
rest are time-driven (TD). DSMC particles outside the 
marked region are treated more efficiently, using time- 
stepping and without any neighbor searches, while the 
DSMC and MD particles inside the marked region are 
treated as in MD (with the addition of stochastic colli- 
sions among DSMC particles). In cases when the solute 
particles are much larger than the solvent particles, NNLs 
are used with a special technique called bounded sphere- 
complexes [5] to handle neighbor searches for the large 
particles. Whenever DSMC particles transfer from a TD 
to an ED cell they are inserted into the event queue with 
an immediate update event, and if needed their neigh- 
borhood region and NNL is constructed. Similarly, when 
DSMC particles transfers from an ED to a TD cell they 
are removed from the event queue and their neighbor- 
hood region and NNL is destroyed. 

The time-driven particles are updated together with a 
fixed time step At. These time step events are treated 
as a special external event and thus processed in correct 
time order with the events scheduled for the particles in 
the combined MD/DSMC region. Stochastic collisions 
are only processed at time step events, exactly as in tra- 
ditional DSMC algorithms. 



4- Open Boundary Conditions 

In simulations of polymer chains in solution in three 
dimensions, a very large number of solvent particles is 
required to fill the simulation domain. The majority of 
these particles are far from the polymer chain and they 
are unlikely to significantly impact or be impacted by the 
motion of the polymer chain. These particles do not need 
to be simulated explicitly, rather, we can think of the 
polymer chain and the surrounding DSMC fluid as be- 
ing embedded in an infinite reservoir of DSMC particles 
which enter and leave the simulation domain following 
the appropriate distributions. 



Such open (Grand Canonical) boundary conditions 
(BCs) are often used in multi-scale (coupled) simulations. 
The combination of a partially time-driven algorithm and 
an unstructured (ideal gas) DSMC fluid makes it very 
easy to implement open BCs by inserting DSMC parti- 
cles in the cells surrounding the simulation domain only 
at time-step events, based on very simple distributions. 
At the beginning of a time step event, after possibly re- 
building the cell masks, the time-driven DSMC particles 
are propagated as usual. If there are external cells, (trial) 
reservoir particles are then inserted into the cells that 
are both external and boundary. The trial particles are 
thought to be at a time t — At, and are propagated by a 
time step At to the current simulation time. Only those 
particles that move into a non-external cell are accepted 
and converted into real particles. Following the insertion 
of reservoir particles stochastic collisions are processed in 
each cell as usual. 

Figure [l] provides an illustration of the division (mask- 
ing) of the cells for the simulation of a tethered polymer 
in two and three dimensions. The division of the cells 
into event-driven, interior, boundary and external cells is 
rebuilt periodically during the simulation. This rebuild- 
ing may only happen at the beginning of time steps, and 
requires a synchronization of all of the particles to the 
current simulation time, a complete rebuilding of the cell 
bitmasks, and finally, a re-initialization of the event pro- 
cessing. Importantly, particles that are in purely external 
cells are removed from the simulation and those that are 
in event-driven cells are re-inserted into the event queue 
scheduled for an immediate update event. 



C. First-Passage Kinetic Monte Carlo (FPKMC) 

An exact AED algorithm for kinetic Monte Carlo 
(KMC) simulation of a collection of diffusing particles 
was recently proposed in Ref. [7J. The main difference 
with MD is that the equation of motion of the particles 
is not a deterministic but rather a stochastic ordinary 
differential equations (ODE). Additionally, reactions be- 
tween particles lead to new types of events such as par- 
ticles changing species, appearing or disappearing. This 
complicates the implementation but is conceptually sim- 
ple to incorporate into the basic AED algorithmic frame- 
work. 

Assuming that all of the required propagators can be 
obtained in closed form and evaluated to within numer- 
ical precision, the FPKMC algorithm is exact in the 
sense that it simulates a trajectory that is sampled from 
the correct probability distribution. For pure diffusion 
with transport coefficient D, the probability c(Ar, At) 
of finding an isolated point particle at a displacement 
Ar at time At is the Green's function for the (time- 
dependent) diffusion equation, dtc = DV 2 c, with no ad- 
ditional boundary conditions. In three dimensions 

-i/2 r Ar 2 



c(Ar, At) = (in At) v exp 



AD At 
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Figure 1: The partitioning of the domain into interior (I) [either event-driven (ED) or time-driven (TD)], boundary (B), and 
external (E) cells in two (left) and three (right) dimensions for a polymer chain of 25 beads tethered to a hard wall. The cells 
are shaded in different shades of gray and labeled in the two-dimensional illustration. The DSMC particles are also shown. 



Particles that have a finite extent, such as spheres and 
cubes, are easily handled by considering their centroids as 
the diffusing point particles. Particles with oricntational 
degrees of freedom are more difficult to handle and for 
now we focus on the sphere case. 

Assume that the protective region Vi is disjoint from 
all other protective regions and the core Ci is restricted to 
remain within Vi. For point particles, one can show that 
the probability c(Ar, At) , conditional on the fact that the 
particle never leaves the interior of Vi, is again a Green's 
function for the diffusion equation but with the addi- 
tional boundary condition that c vanish on the boundary 
of Vi- A single-particle propagator consists of sampling 
from such a probability distribution c\. For simple pro- 
tective regions such as cubes or spheres relatively simple 
closed- form solutions for c\ exist. The probability dis- 
tribution ci is only valid under the assumption that a 
given particle i remains inside V%. From c\ (specifically, 
the flux Vci on the boundary of Vi) one can also find 
the probability distribution that a particle first leaves Vi 
for the first time at a time t and at position f, i.e., the 
first-passage probability Ji(t,f). This distribution can 
be used to sample a time at which particle i is propa- 
gated to the surface of Vi, and then Vi is updated. 

The basic idea of the First-Passage Kinetic Monte 
Carlo (FPKMC) algorithm is to protect the particles 
with disjoint protective regions (an unprotected parti- 
cle has Vi = Ci) and then use single-particle propagators 
to evolve the system. Typically the protective regions 
would have the same position and shape as the particle 
itself but be enlarged by a certain scaling factor [i v > 1 . 
At some point in time, however, two particles i and j 
will collide and thus cannot be protected with disjoint 
regions. Such nearly-colliding pairs are protected by a 
pair protection region Vij (e.g., two intersecting spheres, 
each centered around one of the particles) . A pair prop- 



agator c2(Arj, Arj, At) is used to either find the first- 
passage time, that is, the time when one of the parti- 
cles leaves Vij, or to propagate the pair conditional on 
the fact that both particles remain inside Vij. Analyti- 
cal solutions can be found by splitting the problem into 
independent diffusion problems for the center of (diffu- 

sional) mass r^ M ^ = (Ar^ + Ai\,)/2 and for the dif- 
ference r[j = Ari — Arj walker (with some additional 
weighting factors for unequal particles). This makes the 
pair propagator ci = c^c^ a combination of two 

independent propagators, and c«j , which, are sim- 
pler to solve for analytically. The difference-propagator 
is typically more difficult to obtain analytically because 
the condition for collision dij = Dij forms an addi- 
tional absorbing boundary for the difference walker, i.e., 
a collision occurs whenever the first-passage propagator 
J2 (fy , t^) samples a point on that boundary. For 
repulsive particles [I] the boundary dij = Dij would be 
reflective (zero-flux) instead of absorbing. 

If the particles are cubes closed-form solutions can eas- 
ily be found for the pair propagators [7], however, in gen- 
eral, two-body propagators are considerably more com- 
plex (and thus costly) than single-body ones. One can 
in fact use approximate numerical propagators in which 
the difference walker takes small hops until it gets ab- 
sorbed at one of the boundaries of its protection. This 
is similar to how one numerically predicts pair collisions 
between non-spherical particles in MD simulations using 
time-stepping for just a single pair of particles [2]. We 
have implemented hopping pair propagators for spheres. 
The only difference with the analytical propagators is 
that the hopping trajectory may need to be reversed later 
if a third particle forces a destruction of the pair protec- 
tion before the scheduled pair event. We have used a 
state-saving mechanism in which the state of the random 
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number generator at the beginning of the pair event is 
saved and later restored if the predicted event is canceled. 
If the predicted event does in fact occur we avoid the cost 
of repeating the same hopping trajectory by saving the 
final state of the walk (i.e., the position of the difference 
walker) when predicting pair events. 



1. AED Implementation of FPKMC 

Geometric near-neighbor searches arc an essential com- 
ponent of the FPKMC algorithm, and the same methods 
(LLCs and NNLs) as in MD are used. Cell transfers are 
not explicitly predicted or processed in this algorithm, 
rather, whenever the position of a particle is updated the 
LLCs need to be updated accordingly. When NNLs are 
used, the collision of Ct with Mi may be sampled exactly, 
or, alternatively, the neighborhood Mi may be updated 
whenever Ci is very close to the inner wall of Mi . We say 
that a particle i is protected against particle j or pair jk if 
Vi is disjoint from Vj or Vjk, similarly for pairs of parti- 
cles. The goal of neighbor searches is to protect a particle 
i against other particles and pairs with the largest possi- 
ble Vi- There is a balance between rebuilding protective 
and neighborhood regions too often and propagating the 
particles over smaller steps, and some experimentation 
is needed to optimize the algorithm and minimize the 
number of neighbor searches that need to be performed. 
Whenever a protection Vi is destroyed, particle i should 
be inserted back into the event queue with t e = t, p = oo, 
v = 0, so that it is protected again right away. 

The main AS steps in Algorithm [T] for FPKMC simu- 
lations are: 

• Steps [2] and [3] may involve the processing of particle 
birth processes, where a particle of a given species is 
introduced into the system to model external fluxes. 
These are typically assumed to occur as a Poisson 
process and therefore the time to the next birth is 
simply an exponentially distributed number, with 
the total birth rate given as the sum of the birth 
rates for each of the species. The birth process may 
be spatially homogeneous or the rate may depend 
on the cell in which the birth occurs. The newborn 
particles are inserted into the queue with p = 0, 
v=-\. 

• Step[6]consists of sampling from c\ or J\ and typ- 
ically also rebuilding Vi as in step[8a| For v > 0, a 
particle decay reaction may be processed, i.e., parti- 
cle i may disappear to produce zero or more "prod- 
uct" particles, which are inserted into the queue 
with p = 0, v = — 1. 

• Step [7] consists of sampling positions and Vj from 
the appropriate distribution: 



— If the event is the decay of i or j, then c. 



(D) 



If the event is the collision of i and j, then 
c (CM) j(d) sampled, and the appro- 
priate reaction (e.g., annihilation, coalescence, 
chemical reaction, etc.) is processed. This 
may destroy i and/or j and/or create new 
particles to be inserted into the queue (with 
p = Q, u=-l). 

If the event is the dissolution of the pair ij, 



then either c. 



(CM) 



and J. 



-1 



ACM) 



are sampled, particle j is inserted back into 
the queue with p = 0, v = 0, and Vi is updated 



and a 



{CM) 



are sampled, and then the decay 



reaction is processed. 



as in step 8a (this may protect the particles i 
and j as a pair again). 

Step [8a] is the primary type of event in step [8] and 
consists of updating Vi- The processing of such 
"virtual" collisions with Vi consists of searching for 
the nearest protection region Vj or Vjk among the 
neighbors of particle i (either using LLCs or NNLs) . 
Particle i is then protected against that nearest 
neighbor. If this makes Vi too small then parti- 
cle j or pair jk is propagated to time t and Vj or 
Vjk destroyed. Finally, Vi, Vij or Vik is constructed 
again, depending on the exact local geometry. 

Step [9] consists of sampling event times from the 
appropriate distributions: 

— In step [9a] J\ is sampled. 

— In step [9b] an exponentially distributed time is 
generated based on the decay rates for species 

Si. 

— In step [9c] J2 is sampled, as well as a decay 
time for each of the two particles, and the ear- 
liest of the three times is selected. 



2. Time-Driven Handling 

Under certain conditions event-driven handling of dif- 
fusion may become inefficient or cumbersome. For exam- 
ple, since all protections are either of single particles or 
pairs of particles, very small protective regions need to be 
used in very dense clusters of particles where there may 
be many nearly colliding triplets of particles. The hops 
taken by the asynchronous algorithm will then become 
just as small as what a time-driven (approximate) algo- 
rithm would use, and it will therefore be more efficient 
to use time-stepping for those particles (and thus avoid 
the cost of event queue operations and also simplify over- 
lap search). Additionally, some particles or boundaries 
(e.g., grain boundaries) may have complex shapes and 
thus their diffusion or interactions with other particles 
may be difficult to treat analytically (even for spherical 
particles it has proven that exact pair propagators are dif- 
ficult to implement). As another example, tightly-bound 
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collections of particles (clusters) may act as a single par- 
ticle that has complex internal structure and dynamics 
(relaxation). 

Just as for the SEDMD algorithm, one can add a time- 
driven component to the asynchronous event-driven FP- 
KMC algorithm. Particles that are time-driven do not 
need to be protected against each other, instead, they 
may be unprotected or they can be protected only against 
event-driven particles. Time-driven particles whose pro- 
tective regions overlap form a cluster and should be 
treated using a synchronous algorithm (this cluster may 
include all time-driven particles). In diffusion problems 
it is often the case that different species have widely dif- 
fering diffusion coefficients and therefore very different 
time-steps will be appropriate for different species. To 
solve this problem, one can use the n-fold (BKL) syn- 
chronous event-driven algorithm inside each cluster. In 
this algorithm, at each event (hop) only particles of a 
single species hop by a small but non-negligible distance 
and the rest remain in place, and more mobile particles 
are hopped more frequently (with the correct relative fre- 
quency) than the less mobile ones. 

Note that a cluster may freely take hops up to the 
time of the next event in the queue without stopping 
and modifying the event queue, since it is known that 
those hops will not be pre-emptied by another event. In 
some situations this may improve efficiency by reducing 
the number of heap operations and also increasing the 
locality of the code by focusing multiple events on the 
same (cached) small group of particles. 



IV. DISCUSSION 

In this section we focus on some of the difficulties in 
deploying AED algorithms in the simulation of realistic 
systems. The best-known difficulty is the parallelization 
of AED algorithms [13 , which we do not discuss here 
due to space limitations. Instead, we will focus on the 
difficulties that make even serial simulations challenging. 
It is important to point out that for problems involving 
hard particles, that is, particles interacting with discon- 
tinuous potentials, event-driven approaches are the only 
exact algorithm. Time-driven algorithms always make 
an error due to the finite size of the time step At, and 
typically At must be much smaller than the actual time 
step between events in order to guarantee that no events 
are missed. 

The most involved aspect of implementing an AED al- 
gorithm for a particular problem is the need for analytic 
solutions for various one- and two-body propagators or 
event predictions. For EDMD, the difficulty is with pre- 
dicting the time of collision of two moving hard particles. 
For hard spheres, this can be done analytically relatively 
easily (but numerical care must be taken [H]). When 
orientational degrees of freedom are involved, however, a 
time stepped ODE-like methodology is needed since an- 
alytical solutions are difficult to obtain [2J. This makes 



collision prediction much more difficult to implement in 
a numerically-stable way and also much more costly. In 
FPKMC, there is a need to analytically construct the 
probability distributions c\, c%, J\ and J2, or at least 
to find a way to efficiently sample from them. These 
distributions are Green's functions for a time-dependent 
diffusion equation inside regions such as spheres, cubes, 
spherical shells, intersection of two spheres, or intersec- 
tion of a cone and a sphere. For time-independent prob- 
lems such solutions can be constructed more easily, but 
for time-dependent problems even the simple diffusion 
equation poses difficulties (analytical solutions are typ- 
ically infinite series of special functions in the Laplace 
domain). Different boundary conditions such as reac- 
tive surfaces require even more analytical solutions and 
tailor-made propagators. The handling of more complex 
equations of motion such as the full Langevin equation 
(which combines convection and diffusion) has not even 
been attempted yet. 

This makes designing a more general-purpose AED 
program virtually impossible. This is to be contrasted 
to, for example, time-driven MD where different interac- 
tion potentials can used with the same time integrator. 
In general, time stepped approaches are the only known 
way to solve problems for which analytical solutions do 
not exist, including two-body problems in the case of 
EDMD or FPKMC. The algorithm used in Ref. [2] to 
predict the time of collision for pairs of hard ellipsoids 
combines a time-driven approach with the event-driven 
one. It does this without trying to combine them in an 
intelligent way to avoid wasted computation (such as re- 
peated trial updates of the position of a given particle 
as each of its neighbors is processed). We believe that 
such an intelligent combination will not only provide a 
more general AED algorithm, but also make the algo- 
rithm more robust numerically. 

More generally, combining event-driven with time- 
driven algorithms is important for efficiency (at a certain 
sacrifice in accuracy) . When the time step is large enough 
so that many events occur within one time step one can 
use the event-driven algorithm in-between the updates in 
the time-driven approach [SJ [T5J [57] . When the time step 
is small, however, events occur sparingly only in some of 
the time steps, and a different methodology is needed. 
We proposed to add a new kind of time step event that 
indicates propagation over a small time step (e.g., for the 
set of particles in pure DSMC cells). The essential ad- 
vantage of event-driven algorithms is that they automat- 
ically adjust to the time scale at hand, that is, that they 
take the appropriate time step without any additional 
input. The real challenge is to use time stepping in an 
event-driven framework in which the time step is adjusted 
accordingly to not waste computation, while still keeping 
the approximations controlled. We have described such 
a framework for the SEDMD and FPKMC algorithms. 
The merging of discrete-event and continuous-time mod- 
els has also been formalized in a rather general simulation 
framework |28j . 
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Another unexplored or barely explored area is that 
of using controlled approximations in AED algorithms. 
Approximate event-driven algorithms have been used to 
handle a variety of processes, however, these often use un- 
controlled approximations. The approximate algorithms 
may reproduce the required (macroscopic) physical av- 
erages just as well as the exact algorithm would, how- 
ever, controls are necessary to validate the simulations. 
Examples of approximations that may be useful include 
ignoring unlikely interactions between certain particles, 
approximate solutions instead of exact propagators (such 
as expansions around the mean behavior), etc. 

Almost all of the AED particle algorithms to date have 
focused on single-particle or pair events. This is possible 
to do for hard particles because exact triple collisions are 
extremely unlikely to occur. However, for more realistic 
models, or when approximations are made, events involv- 
ing clusters of particles may need to be considered. For 
example, a cluster of particles may evolve as a strongly- 
coupled (e.g., chemically bonded) unit while interacting 
with other (e.g., freely diffusing) single particles or clus- 
ters. An additional assumption in most AED particle 
algorithms is that events affect only one or two particles, 
so that the event predictions of the majority of particles 
remain valid after processing an event. In some situa- 
tions, however, there may be global degrees of freedom 
and associated events that affect all of the particles. For 
example, in MD there may be a macroscopic strain rate 
that affects all of the particles, since all of the event pre- 
dictions are invalidated when the strain rate changes. In 
principle the strain rate is coupled back to each of the 
particles, so that every particle collision also changes the 
strain rate (albeit by a small amount). In time-driven 
MD this is no problem since the evolution of the sys- 
tem is synchronous and the strain rate evolves together 



with the particles, however, in event-driven MD such cou- 
pling between all of the particles makes it impossible to 
schedule events efficiently In this work we restricted our 
attention to problems with only short-range interactions 
between the particles. However, many problems of inter- 
est include long-range (e.g., electrostatic) interactions as 
well, and these interactions effectively couple the motions 
of all of the particles. 

Finally, multi-algorithm and/or multi-scale combina- 
tions including an AED component have not been ex- 
plored to our knowledge. As an example, consider 
the simulation of nano-structures during epitaxial film 
growth [8^. At the smallest scales, time-driven (first- 
principles or classical) MD is needed in order to study 
the attachment, detachment, or hopping of individual 
particles or clusters. Once the rates for these processes 
are known, lattice-based KMC can be used to evolve the 
structure more quickly without simulating the detailed 
(vibrational) motion of each atom. At larger scales, the 
continuum-based FPKMC algorithm we described can be 
used to propagate atoms over large distances in lower- 
density regions (across flat parts of the surface). Fi- 
nally, a time-driven continuum diffusion partial differ- 
ential equation solver can be used to model processes at 
macroscopic length scales. Such ambitious investigations 
are a challenge for the future. 
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